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ABSTRACT 



A method is presented for investigating the periodic sig- 
nal content of time series in which a number of signals 
is present, such as arising from the observation of multi- 
periodic oscillating stars in observational asteroseismol- 
ogy. Standard Fourier analysis tends only to be effective 
in cases when the data are perfectly regularly sampled. 
During normal telescope operation it is often the case that 
there are large, diurnal, gaps in the data, that data are 
missing, or that the data are not regularly sampled at all. 
For this reason it is advantageous to perform the analysis 
as much as possible in the time domain. Furthermore, for 
quantitative analyses of the frequency content and power 
of all real signals, it is of importance to have good esti- 
mates of the errors on these parameters. This is easiest to 
perform if one can use linear combinations of the mea- 
surements. Here such a linear method is described. The 
method is based in part on well-known techniques in ra- 
dio technology used in every FM radio receiver, and in 
part on the SOLA inverse method. 
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1. INTRODUCTION 



Low-amplitude solar and stellar oscillations can be con- 
sidered as a linear superposition of normal modes of res- 
onance. In theory one indexes each by a radial order n 
for the mode, and the degree I, and order m of the spheri- 
cal harmonic functions which describe the azimuthal spa- 
tial behaviour of the normal modes. For solar-like os- 
cillations, excited stochastically by turbulent convection, 
there are usually many oscillation modes present at any 
epoch, with comparable amplitudes. If a time series, in 
which there are a number of periodic signals present, is 
sampled perfectly regularly the problem of determining 
the frequencies can be solved in a straightforward man- 
ner using a Fourier transform and there is a single peak 
for each real frequency present in the data. 

If there are gaps in the data or if the data is irregularly 
sampled there will be more than one peak in the Fourier 
spectrum for each signal frequency in the data. The shape 



of this response for each signal frequency that is present 
in the data is known as the 'window function'. Clearly, 
if there are a large number of real signal frequencies in 
gapped or irregularly sampled data, the Fourier transform 
very quickly becomes very complicated and hence it be- 
comes very difficult to determine which are true frequen- 
cies, and which are aliases generated by the sampling. 

There exist a number of techniques for resolving the 
problem of determining amplitudes, frequencies and 
phases of periodic signals in irregularly sampled, 
gapped data. Some attem pt to 'fill gaps' us i ng au - 
toregressive methods (cf. 'Fahlman & Ulrvch I (|1982|), 
Brown & Christensen-Dalsaaard ( 1990), iFbssatetall 
019991) ). Other methods attempt to do a least-squares 
fitting of sinusoids in the time domain (cf. ScarglS 
(I1982I) . IScargle (1989), Koen (19921)) or to define a 
best-fit spectrum in a least-squares sense, but mollified 
to take into account that this probl em is ill-posed sinc e 
it is a fo rm of an inverse problem dStoica et alJ (l2000h . 
IVio etaL(2nnn)). 

Apart from the problems of sampling, there is also the 
problem of noise in the data. In principle a star can pro- 
duce 'intrinsic noise', through the stellar equivalent of 
granulation. The telescope and instrumentation can also 
introduce noise into time-series which in an ideal case 
would amount to the photon noise due to finite integration 
times. Such noise might be white, i.e. independent of fre- 
quency, but in some cases it could be confined to bands 
of frequencies. It therefore makes some sense to 'chop 
up' the full frequency range available, set by the sam- 
pling and length of the time series, into different regimes. 
Otazu et ak. ( 2004.) refer to this as a 'multiresolution ap- 
proach' . 



2. THE METHOD 



In the approach of lOtazu et al.l (l2004ft the series is con- 
volved with every member of a set of Gaussians, with 
successively increasing widths Aj+i = 2Aj and the dif- 
ferences between these successive convolutions are com- 
puted. If the data were regularly sampled each of these 
differences would correspond in the frequency domain to 
having applied a band-pass filter to the data. If all the 
different multiresolution levels were to be co-added the 



original data would be reproduced. The convolution is 
carried out in the time domain by defining weights Wi 
for each datum Yi measured at time ti and then summing 
WiYi for all i. 

This method of separating frequency space into broad 
bands has the convenience of simplicity since convolving 
with Gaussians in the time domain corresponds to multi- 
plying by Gaussians in the Fourier domain and the filter 
parameters are easy to determine. However, the shape of 
the filter corresponding to the differences between lev- 
els is awkward since it is skewed and it reduces the am- 
plitude of the signal everywhere except at one off-center 
frequency in the band. It is preferable to have a filter that 
is flat over the band, at least in the limit that errors ai are 
equal for all i, with a symmetric smooth drop-off at the 
lower and upper limit of the band. One example of that is 
a filter defined by the following convolution : 
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Band-pass filtering, i.e. multiplying in the Fourier do- 
main by this filter function H{ijj — Wc), centered on a fre- 
quency Wc, corresponds in the time domain by convolving 
with the function : 
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h[t) = exp 

nt 



cos U)ct 



(2) 



One can create successive bands that are orthogonal by 
choosing for instance : 

= r 1 fixed 

Ao, 

A^j+i = ajA^j 

= 7,+i = l + ^ with 71 = 0(3) 

The weights Wm,ij are therefore in this case set as : 

A, 
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2.1. local oscillator step 

In order to obtain the frequency content of the time se- 
ries, without performing Fourier transforms on unevenly 
sampled data, one can employ what is known as a Hilbert 
transform. In radio receiver technology this would be re- 
ferred to as mixing the signal with a local oscillator, and 
then applying a low-pass filter. The frequency / of the lo- 
cal oscillator is 'tunable' : the width of the band-pass set 
by the previous step is covered, sampling regularly in fre- 
quency with some sampling step A/. This sampling dis- 
tance A/ is chosen in combination with the width of the 
low-pass filter subsequently appUed, so that some over- 
sampling is performed. 

The low pass filter to apply to the 'mixed signal' is cho- 
sen to be a Gaussian. This is a convenient choice because 
multiplying by a (narrow) Gaussian filter in the frequency 
domain, corresponds to convolving with a (wide) Gaus- 
sian in the time domain. The width A^pf in the time 
domain of what corresponds to a low-pass filter, is lin- 
early related to the parameter A of the multire solution 
level within which the signal is to be frequency analysed. 

One proceeds as follows. For each frequency fk define 
sets of 'local oscillator' weights qi^k, Pi,k '■ 



qi^k = 2cos(27r/fcii) 
Pi,k = -2sin(27r/feii) 

and also define a set of low-pass filter weights Zi^i : 



Zi.l = 
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(8) 



The central times Ti are spaced by some factor of or- 
der unity times Alpf- A^p^is large in order to only 
let through signal in a narrow band, and therefore the 
number of times T; is small : typically of the order of 
the number of nights of observations. One can define a 
cosine-weighted average (-Rcos) and a sine-weighted av- 
erage (i?sin) as follows : 
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where Wj is again a normalisation factor and the x. 
and hj are : 



The three indices for these two averaged quantities refer 
to : 
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The band-pass filtered data R are now calculated as : 

Rij = RtiJ = ''^WTn,ijYTn (6) 



1. j the 'resolution level' of the time series, which cor- 
responds to a broad range in frequencies between 
a lower and upper bound set by a smoothing width 
A,. 

2. k the central frequency fk of the tunable narrow- 
band filter runs from the lower to the upper limit of 
the band j to explore the amplitude of signal in the 
time series at/around each fk. 



3. I the central time on which a broad Gaussian is cen- 
tered, providing a narrow band fiher By taking a 
finite value for the width of this Gaussian (rather 
than infinite) and stepping through the time series, 
a time-frequency analysis is done, since it provides 
a mechanism to consider the frequency content per 
night, or per few successive nights. 

For data covering only a few nights with a large gap dur- 
ing the day-time, the window function has quite strong 
daily sidelobes which complicates analysis of a spectrum 
if there are a number of peaks present, as is common in 
multi-periodic variable stars. It is therefore useful to con- 
sider improvements on the method with the aim of reduc- 
ing sidelobe structure. This may be achieved by choosing 
the low-pass filtering weights Zi^i in a different manner 
than described by (|8|i. 



2.2. Strategy 

Consider again the two summations of Eq. (|9} which are 
defined to be functions of frequency /. The weights are 
now Ci and in principle one wishes to be free to choose 
weights differently depending on the frequency fk of the 
local oscillator. The time 7] on which the low-pass filter 
is centered in Eq. ^ now enters the discussion by ex- 
plicitly setting the reference time for the phase (p, which 
means replacing with ti — Ti. With the weights ( yet to 
be determined two functions of / are defined as follows ; 

Y,C^,kl cos {2Trif-fk)iU~Ti)) = 

i 

Y,C^.klS^^iM.f - fk)iU-Ti)) = X{f)(lO) 

i 

The method has three separate aims to achieve with its 
choice for the weights. These aims may well not be per- 
fectly compatible and therefore require a form of com- 
promise : 



• The first function needs to be a function 

peaked at / = and smoothly dropping away, with 
as httle side-lobe structure as possible. For instance : 
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with as small a width A/ as achievable. 

• The second function X{f) needs to be as close to 
as possible everywhere 

• The errors in the resulting amplitudes need to be 
kept as small as possible. 

In the framework of the S OLA m ethod ^ 
Piicers & Thomcson (1921, iPiiners & Thomosc 
( 119941) ^ the best compromise between these three aims 
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Figure 1. A histogram of the distribution of intervals be- 
tween successive sampling times for an artificial time se- 
ries, divided by the smallest of these intervals. The ver- 
tical dashed line indicates the mean time interval for the 
entire set. The inset shows the part of the distribution at 
small intervals which is off-scale for the larger figure. 



can be achieved by minimising for ^ the quantity 
>V(C) defined by : 

/high 

W(C)^(1-m) / d/[S(/)-T(/)]^ 

/low 

/high 

j d/[X(/)]VA^G,fe0.fcSy (12) 

/low ''^ 

while taking into account the normaUsation constraint : 



/high 



d/cos(27r(/-/fc)(t,-Ti)) 
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= 1 



(13) 

Here the target function T(/) is taken to be a Gaussian as 
in ( II 1> . The 'target function' for X{f) is 0. The matrix 
Yiij is the variance-covariance matrix of the data errors. 
It is usual to assume that the errors on the data are Gaus- 
sian distributed and independent, so that E is a diagonal 
matrix, but this is not essential for the algorithm. The 
methods described in Koen ( 1999) to estimate or model 
error properties from residuals using ARIMA models can 
be applied here as well. The parameters ^ G [0, 1] and 
A G [0, cxd) weight the relative importance of the three 
terms. Their value needs to be set by trial and error. A 
very low choice of A normally leads to unacceptably high 
errors on the inferred A^, whereas a very high value gen- 
erally produces the undesirable effect of increasing width 
or structure in S. Similar arguments hold for /i. 

A point that is crucial for this version of the SOLA algo- 
rithm, is that /low and /high must be finite. This is nec- 
essary because the integrations in Eq. (I12> are over prod- 
ucts of cosines, which causes problems if they are to be 
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Figure 2. The response to a monochromatic signal of 2.5 mHz sampled with the sampling shown in Fig. 1, which is 
the equivalent of the window function for this method. The left and middle columns show the functions E! and X, with 
in addition the target functions (dashed linmes). The right hand column shows the sum of the square of these which 
is the equivalent of the window function for power. The rows are in order of increasing spectral resolution with the 
A/ = 8, 4, 2, 0.69 /LtHz, from top to bottom. The lower and upper frequency limit of the band are set to /low = 2.5 mHz 
and /high = 6 mHz. 
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Figure 3. The window function obtained using the algo- 
rithm of section 2.1 with parameter choices correspond- 
ing to the bottom right panel of Fig. 2. 

integrated over an infinite range. The algorithm is there- 
fore limited in application to bandwidth-limited data. In 
practice this is not a severe problem, since by construc- 
tion of the multiresolution steps, bandwidth limitation is 
achieved. Also, given a finite minimum interval between 
measurements it is known that signals with frequencies 
above the Nyquist frequency cx 1/ min(ti+i — ti) cannot 
be resolved by such data, which also sets an upper limit 
to the frequency. 



3. RESULTS FOR THE SAMPLING OF AN AS- 
TEROSEISMIC CAMPAIGN 



To show how the methods work, the sampling times 
are used of a time series obtained during an asteroseis- 
mic campaign : the sampling times are nearly regularly 
spaced during blocks of time, with large gaps in between 
the blocks. This corresponds to night-time observing for 
a few successive nights, with day-time interruptions as 
well as loss of data due to adverse weather conditions as 
well as some technical glitches. The total number of sam- 
pling times is 1962, the median sampling rate is about 
^ s whereas the mean rate is x~ 236 s. If the se- 
ries had been regularly sampled at the median rate the 
Nyquist frequency would have been ^7.14 mHz. Fig. ^ 
shows the distribution of sampling time intervals between 
successive measurements, clearly demonstrating that the 
time series is not fully regularly sampled. Using these 
sampling times the linear coefficients C, are determined 
and then used on an artificial signal with a frequency 
close to 2.5 mHz to obtain a window function, shown 
in Fig. |2] Also shown for comparison in Fig. |3lis the 
window function obtained from the method of Sect. 2.1. 
The rows in Fig. |2] are in order of increasing spectral 
resolution with the A/ = 8, 4, 2, 0.69 /iHz, from top 
to bottom, where 0.69/iHz corresponds to the resolution 
limit set by the full length of the time series if it were 
regularly sampled over this period. 

Perhaps disappointingly the spectral resolution that can 
be attained over the band is quite poor However, the 
quality of the window function depends on the ratio of 
the resolution over the width of the band : A//A;^. By 
bandwidth filtering of the data in the same way as dis- 



cussed in section 2, but over a narrower band, a better 
resolution can be achieved within the band. In Fig. |4]the 
resolution Ay is set to 2 /zHz and the limits of the band 
are /low — 2.456 mHz and /high = 2.544 mHz i.e. a 
factor of 40 narrower than before. The target sampling 
frequency in this case is 2.49 mHz. Now the FWHM of 
the window function in power is 3.3 /iHz (top right panel 
of Figllll. Even at a target width A/ = 0.69 /iHz the side- 
lobe structure in the window function for power (bottom 
right panel of Fig 13 is much reduced compared to the 
same resolution for the previous case (the bottom right 
panel of Fig.|2jl. 



4. DISCUSSION 



The method described in this paper, for the analysis of 
multiperiodic, irregularly sampled time series, has sev- 
eral potential uses. The first point to note is that the 
method is a time-frequency method formulated in the 
time domain. This means that the frequencies of the sig- 
nals in an irregularly sampled time series can be mon- 
itored as a function of epoch. This allows for the pos- 
sibility to detect changes in oscillation frequencies with 
epoch. Such changes could arise in asteroseismology 
through similar mechanisms that produce changes in os- 
cillation frequencies of the Sun over the solar cycle and 
are therefore of interest to determine. 

The second point to note is that the method is linear The 
algorithm produces a set of linear coefficients on the basis 
of the sampling times alone. These coefficients can then 
trivially be combined with data, and with measurement 
errors, in order to provide the power or amplitude of the 
time series and also an error estimate for this power, and 
even a phase if required. Evidently this linearity is con- 
venient for standard asteroseismic data acquisition which 
provides photometry or velocity. However, it is arguably 
even more useful if the data obtained is more complex 
such as a time series of images. A particular applica- 
tion lies in addressing another fundamental problem in 
asteroseismology, which is the identification of the sur- 
face node line pattern associated with any particular fre- 
quency : the mode identification problem. For stars with 
a single dominant mode, long baseline optical interfer- 
ometry has been suggested as a means to image the flux 
variation over the surface, thus allowing the identification 
of the mode (cf. ijankov et a l. (2001) for a description of 
the relevant technique). In particular the closure phase of 
interferometric signals between three telescopes is sensi- 
tive to symmetries in the flux distribution over the sur- 
face of stars (cf. Domiciano de Souza et al. (2005)). For 
multiperiodic stars such a technique would suffer since 
the superposition of patterns would reduce the closure 
phase response. However, the acquisition of interfero- 
metric fringes is already done repeatedly as a function of 
time for operational reasons. By combining closure phase 
signals at the appropriate phase in time one could select 
a single frequency to image, thus restoring the closure 
phase signal strength. This technique might be referred 
to as 'stroboscopic interferometry' . The closure phase 
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Figure 4. The effect of reducing the width of the broad band over which the frequencies are sampled on the window 
function. Here the lower and upper frequency limit of the band are set to /low — 2.456 mHz and /high = 2.544 mHz. 
The target sampling frequency is 2.49 niHz. 



is a somewhat more complex form of data than time se- 
ries of velocities or fluxes, and carrying out stroboscopic 
interferometry in practice is facilitated by being able to 
carry out the time-filtering using the coefficients gener- 
ated with this SOLA based algorithm. Appropriate tests 
of the algorithm for this application are in progress (Pij- 
pers, in preparation). 

A third point to note is that this method can also be 
adapted to streamline procedures for model fitting in as- 
teroseismology. The Gaussian form of the target function 
chosen here is intended to obtain measures of the oscilla- 
tion power or amplitude as a function of frequency. How- 
ever, in fitting models to a given star it is not just a single 
frequency but the entire pattern of frequencies for which a 
best match needs to be found. This can be achieved more 
directly by calculating the expected oscillation spectrum 
for each model and taking for the target function T{f) 
not a single Gaussian, but a sum of Gaussians centered 
on model frequencies /„;, in which the nl are the radial 
mode order n and azimuthal degree I, instead of the sin- 
gle /fe. The coefficients obtained in this way, when com- 
bined with the data, will produce a maximal response for 
that model which has a set of frequencies that is the clos- 
est match to the data. Further tests of this are necessary to 
establish the error propagation properties of such an ap- 
proach, including the influence of systematic effects (Pij- 
pers, in preparation). 
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